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ABSTRACT 

The problem of detecting Sunyaev-Zel'dovich (SZ) clusters in multifrequency CMB ob- 
servations is investigated using a number of filtering techniques. A multifilter approach 
is introduced, which optimizes the detection of SZ clusters on microwave maps. An 
alternative method is also investigated, in which maps at different frequencies are com- 
bined in an optimal manner so that existing filtering techniques can be applied to the 
single combined map. The SZ profiles are approximated by the circularly-symmetric 
template t{x) = [1 -I- (x/r^)^]^'^, with A ~ ^ and x = \x\, where the core radius 
Tc and the overall amplitude of the effect are not fixed a priori, but are determined 
from the data. The background emission is modelled by a homogeneous and isotropic 
random field, characterized by a cross-power spectrum Pviviil) with q=\q\. The fil- 
tering methods are illustrated by application to simulated Planck observations of a 
12.8° X 12.8° patch of sky in 10 frequency channels. Our simulations suggest that the 
Planck instrument should detect ~ 10000 SZ clusters in 2/3 of the sky. Moreover, we 
find the catalogue to be complete for fiuxes S > 170 mJy at 300 GHz. 

Key words: methods: analytical - methods: data analysis - techniques: image pro- 
cessing - cosmology: cosmic microwave background - galaxies: clusters 
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1 INTRODUCTION 

The detection and characterisation of the Sunyaev- 
Zel'dovich (SZ) effect is currently an area of considerable 
interest in millimetre and sub-millimetre astronomy. The 
SZ effect distorts primordial cosmic microwave background 
(CMB) radiation in the direction of galaxy clusters due to 
inverse Compton scattering of CMB photons by the hot in- 
tracluster plasma. In the context of the ESA Planck mis- 
sion (and other CMB experiements with high sensitivity and 
high angular resolution), one must correct for this distortion 
(and for other foreground contaminants) in order to study 
the primordial anisotropics in the CMB. Nevertheless, the 
SZ effect itself is of considerable cosmological interest. The 
effect in individual clusters can be used to study the intr- 
acluster medium. Moreover, since the amplitude of the ef- 
fect does not depend on redshift, it can also be used as a 
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means for detecting clusters that would otherwise be unob- 
servable. Perhaps most importantly for cosmology, a cata- 
logue of the SZ effects in a large number of clusters, could 
be used to place constraints on cosmological parameters in- 
dependently of those resulting from analysis of primordial 
CMB anisotropics. Indeed, a goal of the forthcoming Planck 
mission will be to produce a full-sky catalogue of SZ effects 
containing several tens of thousands of galaxy clusters. It is 
therefore crucial to have a robust and efficient method for 
detecting and extracting the SZ effect from multifrequency 
maps of the CMB. 

The issue of component separation using multifre- 
quency CMB maps has been thoroughly studied in the lit- 
erature. Proposed methods include Wiener filtering (WF, 
Tegmark & Efstathiou 1996, Bouchet et al. 1997), the 
maximum-entropy method (MEM, Hobson et al. 1998, 
1999), fast independent component analysis (FastICA, 
Maino et al. 2001), Mexican Hat Wavelet analysis (MHW, 
Cayon et al. 2000, Vielva el al. 2001a), matched filter anal- 
ysis (MF, Tegmark & Oliveira-Costa 1998) and adaptive 
filter techniques (AF, Sanz et al. 2001, Herranz et al. 2001a, 
Herranz et al. 2001b). A non-parametric Bayesian approach 
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to detecting SZ clusters has also recently been proposed by 
Diego et al. (2001b). A comparison between these methods 
is difficult because they have different specific purposes, use 
different sets of assumptions and the quality of the separa- 
tion varies under different circumstances. For example, the 
MHW is designed to detect compact sources with a Gaus- 
sian profile (such as point sources convolved with a Gaussian 
beam) whereas other methods are better suited to deal with 
diffuse components such as dust and synchrotron emission. 
This suggests the possibility of combining several of these 
methods in order to improve the component separation, for 
example MEM + MHW (Vielva et al. 2001b). 

In general, a component separation method that uses 
all the available information will be more powerful than one 
that does not assume any prior knowledge about the data. 
The MEM, for example, produces excellent results when the 
power spectra and the frequency dependences of the com- 
ponents are well known. Unfortunately, if these assumptions 
about the data are incorrect, errors may arise that would 
affect the separation of several (or all) components. This 
is particularly dangerous in methods that perform the sep- 
aration of all the components simultaneously (WF, MEM, 
FastlGA). The opposite approach is to use a robust method 
that makes as few assumptions as possible about the data 
and aims to separate out just a single physical component. 
An example of the latter approach is the non-parametric SZ 
detection method given by Diego et al. (2001b), in which 
only the well-known frequency dependences of the SZ ef- 
fect and the CMB are assumed. In general, however, such 
methods are less powerful than ones which assume a greater 
degree of (correct) prior knowledge about the data and the 
physical component of interest. Clearly, some compromise 
between robustness and effectiveness must be made when 
proposing a component separation method. 

Filtering techniques (such as MHW, MF and AF) are 
single component separation methods that use some of the 
available information (i.e. the spatial structure of the com- 
ponent to be detected and the power spectrum of the combi- 
nation of the other 'background' components), and generally 
reduce inaccuracies in the separation that may appear due 
to error propagation in methods that separate all compo- 
nents simultaneously. The MHW assumes a specific shape 
(a Gaussian) for the component to be separated ('sources') 
and, given the power spectrum of the background (that can 
be directly estimated from the data) , finds the optimal scale 
of the filter. The MF generalises by allowing more general 
spatial profiles for the component (usually discrete objects) 
to be separated. Adaptive filters put additional emphasis on 
the characteristic scale of the sources in order to reduce fur- 
ther the number of spurious detections. Although one usu- 
ally assumes spherical symmetry of the sources to be de- 
tected, it is not a general requirement of filtering method 
and filters can be easily generalised to detect asymmetric 
features. 

In this paper, we discuss the generalisation of filtering 
techniques, in particular MF and AF, to the case of multi- 
ple data maps corresponding to different frequency channels. 
Multifrequency information can be used both to increase the 
signal of the sources and to reduce the contribution of the 
background ('noise'). This information can be used prior 
to the filtering of the images (by using the correlations be- 
tween the different channels to find an optimal combination 



of channels that maximises the signal-to-noise ratio of SZ 
clusters) or can be included directly in the construction of 
the filters. The structure of this paper is as follows. In sec- 
tion ^ we describe the formal aspects of the methods that 
make use of multifrequency information proposed in this pa- 
per. In section ^we describe the simulations we made to test 
the multifilters. Section ^ is devoted to summarise and dis- 
cuss our results. 



2 MULTIFREQUENCY FILTERING 

There are two different approaches we can follow to include 
the frequency dependence of a signal in the filtering of multi- 
channel data. On the one hand, we can filter each chan- 
nel separately, but carefully taking into account the cross- 
correlation between the different channels and the frequency 
dependence of the signal in order to obtain an output set of 
filtered maps that, when added to each other, is optimal for 
the detection. This philosophy gives birth to the multifilter 
method. On the other hand, one can use the information 
about correlations and frequency dependence before filter- 
ing in order to find the optimal combination of channels that 
maximizes the signal-to-noise ratio of the sources, and then 
use a filter on the optimally combined map. This leads to the 
design of a combination method plus a single filtering. These 
two methods are discussed in detail below. We first, however, 
define our model of the multifrequency observations. 

2.1 Model of the observations 

Let us consider a set of 2-dimensional images (maps) with 
data values defined by 

y^{x) = f^s„{x) + n„{x), x = \x\, u=l,...,N, (1) 

where A'' is the number of maps (or number of frequencies) . 
The first term on the RHS represents the signal in each 
frequency map due to the thermal SZ effect in clusters, and 
the generalised noise nv{x) corresponds to the sum of the 
other emission components in the map. 

For a map of any reasonable size, we would expect the 
presence of several SZ clusters. To illustrate more clearly 
the construction of multifilters, however, we assume that 
the signal is due to a single SZ cluster at the origin £ = of 
the map (the generalisation to several clusters distributed 
across the map is straightforward). In particular, we as- 
sume a spherically-symmetric /3-profile for the cluster elec- 
tron number density 

ne(r) oc [l^{x/r,f]-i^, 

where Tc is the core radius and we adopt the standard value 
13 = 2/3. One trivially obtains that the two-dimensional mi- 
crowave decrement from such a cluster has the form At(x), 
where A is the amplitude of the effect and t[x) is the spatial 
template 

^(^) = [l + (^)2]l/2- (2) 

At each observing frequency this template is convolved with 
the corresponding antenna beam, which we assume to be a 2- 
D Gaussian of dispersion 6^ , to produce the convolved profile 
Tvix). Thus, in (^, Si,{x) = Ar^ix). The quantity in (^) is 
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the frequency dependence of the SZ effect, normahsed such 
that = 1 at the fiducial frequency Vf. Ifence A is the 
true ampUtude of the SZ effect at the frequency Vf. 

The background n^{x) is modelled as a homogeneous 
and isotropic random field with average value {n^{x)) — 
and cross-power spectrum P^iv^io) il = \<l\) defined by 

where n^{q) is the Fourier transform of nu{x) and 5^ is the 
2-D Dirac distribution. 



2.2 The mutlifilter method 

We now consider the first of two approaches to the filter- 
ing of mutlifrequency data, namely the multifilter. In fact, 
multifilters themselves can be constructed in several differ- 
ent ways. We consider here the two main possibilities, which 
are the scale-adaptive multifilter and the matched multifilter. 



2.2.1 Scale-adaptive multifilter (SAME) 

The idea of an optimal scale-adaptive filter was recently pro- 
posed by Sanz et al. (2001), which we now generalise to 
multifrequency data. One begins by introducing, for each 
observing frequency, a circularly-symmetric function '^^^{x). 
From this function, one can define a family of filter functions 



(3) 



where the 3 parameters Ru and h define a scaling and a 
translation respectively. For any particular values of these 
parameters, we define the filtered map at the i^th observing 
frequency by 



(4) 



'Wv{Rv,h)= J dxy,j{x)'^v(x\ Rv,b). 

and the 'total' filtered map as 
■w{Ri, ...,Rn,b) = 'y^ Wi,{R^,b). 



(5) 



The convolution ^ can be written as a product in 
Fourier space, in the form 



Wi,{R^,b) = I dqe ^'''^ y^{q)^^{R^q), 



(6) 



where yu{q) and ipv{q) are the Fourier transforms of yv{x) 
and Tpuix), respectively. A simple calculation then shows 
that the expectation value of the filtered field at the ori- 
gin & = 0, is given by 



{w^{Rv,Q)) =2TiAf^ I dqqT^{q)xlj^{R^q) 



(7) 



where the ensemble-average is over realisations of the back- 
ground emission (x) and the limits in the integrals go from 
to (xj. Similarly, one finds that the variance of the total 
filter field (j^ is given by 

al,{Ri, R„) = {w^{Ri, Rn, b)) - {w{Ri,..., R„, b)f 



= 2-K 



dq q ^2 (g) -iAfi (-R-^i <?) V'i'2 <?) ■ (8) 



We choose the filter functions '4>v{q) to optimise the 
detection of the cluster at the origin, taking into account 
that the source has a bell shape with a single characteristic 
scale Rov in each map. We define the optimal scale-adaptive 
multifilter (SAMF) to be that which ensures the following 
conditions are satisfied: 

(i) w{Roi, Ron, 0) is an unbiased estimator of the am- 
plitude of the source, so {w{Roi, Ron, 0)) = A; 

(ii) the variance of w{Ri , -Rn, 6) has a minimum at the 
scales Roi, Ron, i.e. it is an ejfictent estimator; 

(iii) Wv{Rv,b) has a maximum at (7?oi/,0). 

The multifilter satisfying these conditions is given by the 
matrix equation 

^{q)=p-\aF + G), (9) 

where we have introduced the column vectors '4>{q) = 
[^i,{R^q)\, F = [fvT^], and G = \p,^l3v], where 



d In Tv 
din g 

Also, is the inverse of the matrix P = [Pi^iv2i'l)]t ^^'^ 
the quantities a and are given by the components 

a = (A-')oo, P. = {A-^).o, (10) 
where A is the {1 -\- n) x (1 + n) matrix with elements 

Aoo = J dqF'P'^F, Ao^ = J dq^i,{F'p-')^, (11) 
A,o= / dq^i,{p-^F),, A,,,= / dqfi,fi,,{p-^)^,,. (12) 



The quantity ip{q) in (^) is called scale-adaptive multifil- 
ter, extending the concept considered in our previous paper 
(Sanz et al. 2001) for a single map. Also, from (^, we see 
that the variance can be written as 



dqi>^Pi> 



(^oo)-' 



l + 2y dqF^p-^G -h J dqG^p-'^G. 



An interesting special case is when P is a diago- 
nal matrix, i.e. there is no cross-correlation between the 
backgrounds in the different frequency maps, so P^^i = 
5vv'Pv[l)- 111 I'^i'^ case, the multifilter is given by 



where 



6. = 



2ii 



dq 



P ' 



fl 



dq — [2 -I- 

Pu \ dlnq 

,^T^ /„ , dlnr^ 

di? — 2 -f — 

P„ \ dlnq 



(13) 

(14) 
(15) 
(16) 



If, additionally, one assumes that the backgrounds are white 
noise, i. e. the P,^ are constants and the resolution is the same 
in all maps, producing the convolved profile T{rc,8), then 
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where 

N = 



J u -'- 1/ 



H' 



b_ 
H 



dlrir 
ding 



H 



dq T 



din g 

dlnr 
ding 



(17) 

(18) 
(19) 
(20) 



2.2.2 Matched multifilter (MMF) 

If we do not assume condition (iii) in the definition of the 
SAMF in the previous subsection, one obtains another type 
of multifilter after minimization of the variance condition (ii) 
subject to the constraint (i). The multifilter satisfying these 
less restrictive conditions is given by the matrix equation 



(21) 



where F is the column vector F = [/^r^] and is the 
inverse matrix of the cross-spectrum P. This will be called 
matched multifilter extending the concept usually considered 
for a single map. In this case, the variance of the filtered field 
is given by 



dqip^Pip ■ 



(22) 



In the special case where there is no cross-correlation 
between the backgrounds in the different maps, so P^^i = 
5^^iPv{q), the multifilter is given by 



Pu' 



(23) 



If additionally, one assumes that the backgrounds are white 
noise, i.e. the Pu are constants, and the resolution is the 
same in all maps, producing the convolved profile rifcO), 
then 



fu-tpAq) 



J u -'- f 



N : 



(24) 



A similar result has been independently developed by 
Naselsky et al. (2001). 

2.3 The combination method 

Our second approach to the filtering of mutlifrequency data 
is to perform a two stage process in which one first constructs 
some optimal linear combination of the individual frequency 
maps, and then applies a single filter to this combined map. 

The linear combination of the original frequency maps 
is constructed so that it maximally 'boosts' objects with a 
given spatial template and frequency dependence above the 
background. In general, the combination map is given by 



(25) 



where Ci, is some set of weights (derived below). From (|l|) 
we see that this combined map can be written as 



y(x) = At{x) + e{x), 



(26) 



where the 'combined' template t{x) and background e{x) are 
given simply by 



t{x) 



E 



f„T„{x), 



e{x) 



E 



Ci,n„{x). 



(27) 



The optimal values of the weights are found by maximis- 
ing the quotient Q = At{0)/a^, where (jj is the dispersion 
of the combined noise field e{x). Thus, one is maximising 
the 'detection level' of the object of interest above the back- 
ground. 

In practice, it is easier (but equivalent) to maximise Q^. 
It is straightforward to show that, in the numerator of Q^, 

2 



[t(0)Y = c'Gc, 
where 

c = (c. ) , G = (G,,/ ) , = /.T. (0)/.' r,, (0) . 
Similarly, in the denominator of Q^, 
(jj = c Mc, 
where 

c= (c^), M = (M^^,/), M^^, = {nAx)n^,{x)). 



(28) 
(29) 
(30) 
(31) 



Thus one must maximise, with repect to the weights c^, the 
quantity 

,2 c*Gc 



Q 



(32) 



c*Mc 

The solution is easily found to be the eigenvector c associ- 
ated to the largest eigenvalue A of the generalized eigenvalue 
problem 



(G-AM)c = 0. 



(33) 



Once we have constructed the combined template t{x), 
given in (psj), with the optimal weights c given by (^), we 
can apply a either a scale-adaptive filter or matched filter to 
the single combined image given by (pij). The appropriate 
form of the scale-adaptive filter for a single image is given 
by Sanz et al. (2001). The appropriate form of the matched 
filter for a single image is given by 

1 t{q) 



^{q) = 



where 



2na P{q) ' 



a = I dqq 



t\q) 



P{q)=c'Nc, N=[P,,,{q)], 



(34) 



(35) 



where t{q) and P{q) are the Fourier transform of t{x) and 
the power spectrum of the noise e, respectively. 

2.4 Comparison of methods 

We have introduced two different methods for filtering mut- 
lifrequency data, namely multifilters and the combination 
method. Moreover, each method can be implemented for two 
different kinds of filter (scale-adaptive filters and matched 
filters). This makes four different possible ways to filter the 
data, and we would expect each method to have advantages 
and shortcomings with respect to the others. 

From the methodological point of view, the main differ- 
ence between the combination method plus a single filter and 
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Table 1. Technical characteristics of the 10 simulated Planck 
channels. Column two lists the FWHM assuming a Gaussian 
beam. Column three shows the pixel size in arcmin. Column four 
lists the fractional bandwidths of each channel. In column five 
the instrumental noise is AT {^iK) per resolution element for 12 
months of observation is given. 



Frequency 


FWHM 


Pixel size 


Fractional 


^ noise 


(GHz) 


(arcmin) 


(arcmin) 


bandwidth 










(/\u/v) 




30 


33.0 


6.0 


0.20 


AA 


44 


23.0 


6.0 


0.20 


6.5 


70 


14.0 


3.0 


0.20 


9.8 


100 (LFI) 


10.0 


3.0 


0.20 


11.7 


100 (HFI) 


10.7 


3.0 


0.25 


4.6 


143 


8.0 


1.5 


0.25 


5.5 


217 


5.5 


1.5 


0.25 


11.7 


353 


5.0 


1.5 


0.25 


39.3 


545 


5.0 


1.5 


0.25 


400.7 


857 


5.0 


1.5 


0.25 


18182 



the multifilter method is that the first one uses the multi- 
frequency information to give the optimal starting point for 
the filter, while multifilters produce the optimal ending point 
after filtering. In that sense, multifilters are more powerful 
than a single filter applied to an optimally combined map. 
The cost of this higher efficiency, however, is an increase in 
the complexity of the filters and therefore of the computa- 
tional time required to perform the data analysis. For exam- 
ple, in the application to simulated Planck observations in 10 
frequency channels described below, the multifilter method 
requires ~ 10 times more CPU time than the combination 
method (because the last one filters only once). 

Regarding the type of filter, matched filters give higher 
gains while adaptive filters reduce the number of errors in 
the detections (spurious sources). This is due to the fact 
that the optimal scale- adaptive filter includes a constraint 
(numbered (iii) above) about the scale of the sources that is 
not present in matched filter design. This constraint char- 
acterises with more precision the sources but also restricts 
the minimisation of the variance in the filtered maps (that 
is, the final gain). 



3 SIMULATIONS 

In order to test the multifrequency filtering methods out- 
lined above, we now apply them to simulated Planck obser- 
vations. These simulations mimic the main features of the 
future Planck mission, and are realistic in the sense that 
they include the latest available information about different 
physical components of emission (CMB, Galactic compo- 
nents, SZ effect and extra-galactic point sources) and that 
they reproduce the technical specifications of the different 
Planck channels (pixel sizes, antenna beams and noise lev- 
els). Table ^ shows the assumed observational characteristics 
of the simulated maps. The simulations were performed in 
patches of the sky of 12.8° x 12.8°. However, the method 
can be easily extended to the full sphere. 

The CMB simulation was generated using the Cj's pro- 
vided by the CMBFAST code (Seljak & Zaldarriaga, 1996) 



for a spatially-fiat ACDM Universe with Q.m = 0.3 and 
SIa = 0.7 (Gaussian realization). The Galactic emission is 
assumed to originate from four different components: ther- 
mal dust, spinning dust, free-free and synchrotron. Thermal 
dust was simulated using the template given by Finkbeiner 
et al. (1999). This model assumes that dust emission is 
due to two grey-bodies: a hot one with a dust temperature 
ri5°' ~ 1&.2K and an emissivity a''"' ~ 2.70, and a cold one 
with r^°" ~ QAK and a'^"" ~ 1.67. These quantfties are 
mean values. 

For the free-free template we used one correlated with 
the dust emission in the manner proposed by Bouchet et al. 
(1996). The frequency dependence of the free- free emission 
is assumed to vary as 7^ oc z^""'^^, and is normalised to give 
an rms temperature fluctuation of 6.2fiK at 53 GHZ. 

The synchrotron spatial template has been produced 
using the all sky map of Fosalba & Girardino 0, which is an 
extrapolation of the 408 MHz radio map of Haslam et al. 
(1982), from the original 1 deg resolution to a resolution of 

5 arcmin. We further extrapolated the small-scale structure 
to 1.5 arcmin following a power-law power spectrum with an 
exponent of —3. The frequency dependence is assumed to be 
Ii, oc and is normalized to the Haslam 408 MHz map. 
We include in our simulations information on the changes 
of the spectral index as a function of electron density in 
the Galaxy. This template has been made by combining the 
Haslam map with the Reich & Reich (1986) map at 1420 
MHz and the Jonas et al. (1998) map at 2326 MHz, and can 
be found in the FTP site refered to above. 

We have also taken into account the possible Galactic 
emission due to spinning grains of dust, proposed by Draine 

6 Lazarian (1998). This component could be important at 
the lowest frequencies of Planck (30 and 44 GHz) in the 
outskirts of the Galactic plane. 

The extragalactic point source simulations follow the 
same cosmological model as that assumed for the CMB 
simulation and corresponds to the model of Toffolati et al. 
(1998). The thermal SZ effect was made for the same cosmo- 
logical model. The cluster population was modelled follow- 
ing the Press- Schechter formalism (Press & Schechter 1974) 
with a Poissonian distribution in 9 and (j>- The simulated 
cluster population fits well all the available X-ray and opti- 
cal cluster data sets (see Diego et al. 2001a for a discussion). 
The different components used for the simulation are shown 
at 300 GHz in figure ^ Figure ^ shows the simulated chan- 
nels taking into account all the components, as well as the 
antenna beam effect and the instrumental noise. 

The simulation described above is the same used by 
Diego et al. (2001b). We have chosen this particular simu- 
lation in order to compare results with that work under the 
same conditions. 



4 RESULTS AND DISCUSSION 
4.1 Testing the methods 

Before applying the methods presented in section |^ to the 
simulations described in section ^, let us illustrate how they 

"t ftp:/ /astro. estec.esa.nl/pub/synchrotron 
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Figure 1. Components present in the simulation at 300 GHz. From left to right and from top to bottom the components are: CMB, 
kinetic SZ, thermal SZ, Galactic dust, free-free and synchrotron. The spinning dust component has the same spatial template as Galactic 
thermal dust emission. The units of the maps are AT/T. 



work in a simplified case. For the sake of simplicity, the sim- 
ulations include the same foregrounds and technical features 
as the simulations described in section ^ except for the fact 
that the spinning dust component was not included. Spin- 
ning dust is the weakest component and its contribution is 
almost negligible in most Planck channels. The other main 
difference with the realistic simulations is the thermal SZ ef- 
fect itself. Instead of the realistic SZ clusters, we simulated 
200 clusters of the same size {vc = 0.5 pixels), and ampli- 
tudes uniformly distributed between and the maximum 
amplitude of the clusters belonging to the realistic simula- 
tion. The simulated clusters have the frequency dependence 
of the thermal SZ effect and the radial profile 

where N = ryrdij-^ 4- r-c). Here denotes the core radius 
and Tv is a 'cut scale' that can be interpreted as the virial 
radius of the cluster. The profile given by ( |36[ ) is a modified 
multiquadric profile that behaves as a /3 model for a; << r„ 
and decays quickly for x > > r„ . In order to simplify the sim- 
ulations still further, the kinematic SZ effect is not included. 
In the real case, the contribution of the kinematic SZ effect 
is expected to be ~ 30 times lower than the contribution of 
the thermal SZ effect, thus justifying our approximation. 



The simulated maps were filtered using the four possible 
combinations of filters and methods presented in this work. 
The combination method took about 1 minute of process- 
ing in a 700 MHz PC. The multifilter method took about 
10 minutes of processing on the same computer. Process- 
ing times were slightly higher for the case of adaptive filters 
than in the case of matched filters. To compute the filters 
and the combinations of the maps the low resolution chan- 
nels were re-binned to the resolution of the highest resolution 
channel (1.5 arcmin). Figure ^ shows the filters (in Fourier 
space) employed in the multifilter method. The adaptive fil- 
ters are represented by solid lines, whereas matched filters 
are represented by dotted lines. Filters, specially those used 
in low frequency channels, are quite complicated and show 
several peaks at different wave numbers q that correspond 
to the different scales the filters are trying to identify on 
the images. Conversely, the filters used after the combina- 
tion method are rather simple and are shown in figure ^ 
The differences between matched and adaptive filters in the 
combination method are small. This indicates that there will 
be few differences in the results of both filters. In the case 
of the multifilters the differences are more pronounced, in 
particular for the 217 GHz channel. 

The results from the test can be summarised as fol- 
lows. The number of detections is higher in the multifilter 
case. This is not surprising since by definition the multifil- 
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Figure 2. Simulated Planck channels. Each map corresponds to the same 12.8° X 12.8° area of the sky at the frequency of the channel. 
The represented channels are, from left to right and from top to bottom, 30, 30, 44, 70, 100 (LFI), 100 (HFI), 143, 217, 353, 545 and 
857 GHz. The units of the maps are AT/T. 



ter method is more powerful than the combination method. 
In particular, if one raises the detection threshold, the dif- 
ference between the two methods increases. As an example, 
table 1^ shows the results of one of the tests. After filter- 
ing, the sources were detected by looking for peaks above 
a 4(T threshold and then compared with a catalogue of the 
original simulated clusters. By comparison, in the same sim- 
ulations, if the detection threshold is raised to Scr (this case 
is not shown in table |^ the matched multifilters produce 
90 detections whereas the matched filter in the combination 
method gives only 78. 

A detected peak is considered a spurious detection if 
the distance to the closest object in the original catalogue is 
greater than 1.5 pixels. The number of spurious detections is 
strikingly low in all cases. Nevertheless, the adaptive filter 
seems to work better than the matched filter (for exam- 
ple, in table ^ the adaptive filter gives spurious detections 
instead of 1 with the combination method). Of all the com- 
ponents present in the simulations, point sources are the 
most likely to produce spurious detections due to the simi- 
lar scale of sources and SZ clusters. However, the frequency 
dependence of the SZ effect greatly reduces the probability 
of contamination. 

The position of the sources was determined with errors 
below the pixel size (1.5 arcmin). However, the error in the 
determination of the amplitude is quite large (~ 35%). Ta- 
ble 1^ shows that this error is due to a systematic bias. To ex- 
plain this bias let us consider the case of matched filter. The 
filter normalization is given by integral a = dqqr^/P, 
being r the source profile in Fourier space and P the power 
spectrum of the background. However, when we analyse a 
patch of the sky we do not have information about the 



power spectrum at all the wave numbers q. An image di- 
vided into pixels is limited by a minimum value qmin and a 
maximum qmax- Therefore, the normalization we calculate 
is a' — /^^"'"'^ dqqr^/P. For the case of images with the same 
size and pixel scale as our simulations, and a multiquadric 
profile T with Vc = 0.5 pixels, the normalisations can differ 
by 10% for a background spectral index 7 = 1 to 60% for a 
background spectral index 7 = 3. The case for adaptive fil- 
ters is more complicated but qualitatively similar. This bias 
is independent of the source fiux and can be calibrated using 
simulations. 

Figure ^ shows the contribution of each channel to the 
source amplitude estimation for each method. The channels 
with larger contribution are the 143 GHz and the 353 GHz. 
This is not surprising since they are the channels with more 
SZ contribution. The 100 GHz HFI channel has more con- 
tribution than the 100 GHz LFI because of its better signal- 
to-noise ratio. The combination method puts more empha- 
sis in the 100 GHz and less in the 143 GHz channels than 
the multifilter method. Contributions from the 30 GHz, 217 
GHz and 857 GHz channels are negligible. However, that 
does not mean that these channels do not contribute to the 
filter construction. For example, if we repeat the analysis 
with only the 5 channels with more contribution (70, 100 
LFI, 100 HFI, 143 and 353 GHz), the number of detections 
reduces by ~ 10% and the number of spurious detections in- 
creases to 6 (at 4o" detection threshold) for the combination 
method case (for the two filters) and to 9 (adaptive filter) 
and 5 (matched filter) for the multifilter method. Even the 
217 GHz channel is important; repeating the analysis with 
all channels except for the 217 GHz increases the error in 
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Figure 3. Filters used to test tlie multifilter metliod. Adaptive filter {solid line) and matched filter (dotted line) are represented in 
Fourier space for each frequency channel. 



the determination of the ampUtudes and also increases the 
number of spurious detections at low a thresholds. 

We can summarise the conclusions of the test as follows: 

• Multifilter method is more powerful in the detection 
and estimation of cluster parameters. 

• Combination method is faster than the multifilter 
method. 

• Multifrequency information reduces the number of spu- 
rious detections. Therefore, it is not critical to use an adap- 
tive filter to that end. The matched filter allows one to detect 
more sources. 

• Some channels contribute more than others to the anal- 



ysis, but all of them carry valuable information; the analysis 
should include all the available data. 



4.2 Results for realistic simulations 

Taking into account the insights provided by the test pre- 
sented in the last sub-section, we are now prepared to con- 
front the analysis of realistic simulations. These realistic sim- 
ulations have been described in section ^ The main differ- 
ence with respect to the previously performed test is that 
in the realistic case clusters have different sizes that are not 
known a priori. Herranz et al. (2001b) proposed a method 
to accommodate with this problem; 
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Table 2. Results of the test for 10 frequency channels and clusters of equal size {jc = 0.5 
pixels). The detection was performed over the 4(t threshold in all cases. The filtering method is 
listed in the first column. Second column shows the number of true detections found over the 
4(T detection threshold. The third column indicates the number of spurious sources detected 
in each case. The fourth column lists the mean error in the determination of the position 
of the detected sources. Column five lists the mean relative bias in the determination of the 
amplitude, defined as bA = 100 < {A — Aq)/Ao >. Column six shows the mean relative error 
in the determination of the amplitude, defined as Ea = 100 <\ A — Aq \ /Aq >. 
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Figure 4. Filters used to test the combination method. The adap- 
tive filter {solid line) and the matched filter {dotted line) are rep- 
resented in Fourier space. 



• Choose a trial core radius Tc and construct the corre- 
sponding filters. 

• Convolve the data with the filters, varying the scales R,j 
of eq. (^. This is performed by substituting tp{q) by tf}{qx), 
where x = R'^/Ri, is simply a dilation factor. 

• In the c ase of scale-adaptive filters, condition (iii) (see 
section 2.2.1 ) implies that, if the assumed value of corre- 
sponds with the true core radius of the cluster, the coeffi- 
cients will be maximum when a; = 1. If this is not true, the 
trial rc is discarded. 

• In the case of matched filters, condition (iii) no longer 
holds. However, a similar criterion can be used: the best 
performance of the filter will occur when the scale of the 
filter and the scale of the cluster is the same. Therefore, if 
the signal-to-noise ratio of the coefficients is not maximum 
for X — 1, the assumed radius rc is not the optimal choice. If 



Figure 5. Contributions of the different channels to the source 
amplitude estimation. The contributions are normalised so that 
the channel with maximum contribution has contribution 1. Pan- 
els corresponding to the combination (matched) and combination 
(adaptive) are the same, since the weights are assigned before 
filtering. However, they are both shown for comparison purposes. 



this condition is verified for more than one value of rc, the 
one with the highest signal-to-noise ratio is chosen. 

• We repeat the process with as many different values of 
Tc as desired. 

In Herranz et al. (2001b), the above method has been suc- 
cessfully applied to single frequency maps containing sim- 
ulated multiquadric profiles and different kind of back- 
grounds. 

The clusters in the realistic simulations have a profile 



r(r) = 



: tan 



(r/rc)2 



1 + (r/rc)2 



(37) 
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where p = /tc is the ratio between the virial radius and 
the core radius of the cluster. This profile is realistic from the 
point of view of the physics of the Sunyaev-Zel'dovich effect. 
The profiles (||) and (||) are almost identical for r < r„. In 
order to calculate the filters in Fourier space, it is easier to 
work with the modified multiquadric (^6|). Therefore, as a 
first approximation we will assume that the clusters can be 
described by the profile (|3^). As we will see, this is a good 
approximation. 

Following the results of the test in the last subsection, 
we choose a matched multifilter to perform the analysis. Af- 
ter applying our method to the simulations we detect the 
clusters by looking for sets of connected pixels above a cer- 
tain threshold. The maximum of these sets give the position 
and the amplitude of the sources. At the Scr level (regions 
with 5 or more connected pixels) we are able to detect 63 
cluster candidates, of which 62 correspond to real clusters. 
The spurious detection appears in one of the borders of the 
image and therefore can be considered as an edge effect. 
The mean error in the determination of the position of the 
clusters was ~ 1 pixel. 

Using the multi-scale analysis we are able to determine 
the core radii of the clusters with a mean absolute error of 
0.30 pixels. The mean bias in the determination of the core 
radii was —0.15 pixels. Since pixelisation effects are expected 
to corrupt structures whose typical scale is much smaller 
than the pixel size, all clusters with correspondingly small 
core radii will not follow the multi-quadric profile. Most of 
the clusters have very small core radii and therefore can 
be considered as point-like sources. To detect these clus- 
ters, a Gaussian profile (corresponding to the beam at each 
channel) was assumed instead of the multiquadric given by 
(p6|). The separation between clusters that are considered 
as point sources and extended sources was set in — 0.4 
pixel, since if is below this limit the FWHM of the mul- 
tiquadric profile is of the order of a pixel. All the clusters 
that were detected as point sources were considered to have 
the (arbitrary) value rc =0.1 pixel. The bias in the deter- 
mination of the core radius mentioned above is due to the 
asymmetry introduced when assigning Vc = 0.1 pixel to all 
clusters whose core radius is less than 0.4 pixel. Figure ^ 
shows the resulting maps after filtering with different scales 
(point source in the left, rc — 1.0 pixel in the middle and 
rc — 2.5 pixel in the right side). Small clusters stand out 
mainly in the left panel, while large clusters are emphasized 
in the right panel. 

Figure |^ shows the detected amplitudes of the 60 bright- 
est detected clusters vs. the true amplitudes. The brightest 
sources in the graph were recovered with relative errors of 
~ 30%. As might be expected, the error increases as the true 
amplitude decreases. For the weakest clusters, the detected 
amplitude tends to a horizontal line. This bias is due to 
the incorr ect determination of the normalization described 
in section 4.1, as well as two additional factors: first, only 
the faint clusters that are overamplified are able to reach 
the detection threshold; second, the detection is performed 
by looking for peaks in the data and, hence, those clusters 
which are by chance enhanced by residual fiuctuations of the 
background (noise) are more likely to be detected. The bias 
due to the normalization should be the same for all sources 
with a given scale and therefore could be calibrated from 
simulations, whereas the second kind of bias affects mainly 



0.0001 

True Amplit\ide ( AT/T 



Figure 7. Amplitudes of the 60 brightest deteeted clusters vs. 
true amplitudes. The dotted line corresponds to the perfect situ- 
ation True Amplitude = Detected Amplitude. 



the weak clusters (that dominate the number counts in this 
case) and can not be easily calibrated. However, the ampli- 
tude estimation could be improved by refining the detection 
method. 

The catalogue of detected sources was complete above 
0.17 Jy (at 300 GHz). Below this limit the completeness 
gradually decreases. There are a few very faint clusters that 
are detected with this method (fluxes of few tens of mjy), 
but the practical detection limit is 40 — 50 mJy (at 300 
GHz). The 62 detections in our small sky patch of (12.8°)'^ 
corresponds to ~ 10000 detections in 2/3 of the sky (the 
area not dominated by the Galaxy). 

The performance of the other filtering techniques pro- 
posed in this work was compared to the performance of the 
matched multifilter for the realistic Planck simulations. For 
the comparison, only the point-like clusters (that dominate 
the number counts) were co nsid ered. The scale-adaptive 
multifilter, as seen in section 



4.1 



performs worse than the 
matched multifilter because their main advantage, that is 
the removal of spurious sources, has already been accom- 
plished by the multi-frequency information included in the 
filters. Since the clusters one is trying to detect are very 
faint, it is the overall gain which is most important. At 
3(7, the scale- adaptive multifilter detects 40 per cent fewer 
clusters than the matched multifilter. Both matched and 
scale-adaptive filtering of an optimally combined map are 
slightly less sensitive than multifilters, producing, in both 
cases, around 20 per cent fewer detections at 3a than with 
the matched multifilter. In other words, in the case of single 
filters, matched and scale-adaptive perform very similarly 
(as was expected from figure^, where both filters look alike). 

The results obtained with the matched multifilter are 
comparable in number of detections and reliability at the 
3cr-threshold to those obtained by Diego et al. (2001b). To 
perform their method it is necessary first to clean the maps 
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Figure 6. Multifiltered data at several scales. Left panel shows the output after filtering with a matched multifilter that considers that 
clusters are point-like. Center panel shows the output after filtering with a matched multifilter for a modified multiquadric profile with 
rc = f .0 pixel. Right panel shows the output after filtering with a matched multifilter for a modified multiquadric profile with Tc = 2.5 
pixel. Note how diff'erent structures are enhanced at different scales. The units of the maps are AT/T. 



carefully. This cleaning includes point source removal with 
a MHW technique (Vielva el al. 2001a), dust substraction 
using the 857 GHz map and CMB substraction in Fourier 
space up to a certain wavenumber limit using the 217 GHz 
map. In order to compare results, we performed the same 
cleaning method and then applied the matched multifilter. 
At 3(T detection threshold, we obtained a 7 per cent more de- 
tections than before. At 2a, we also detected 7 per cent more 
clusters and the number of spurious detections decreased by 
about 12 per cent. This improvement is small given all the 
cleaning work that is necessary to obtain it. Moreover, it 
indicates that the multifilter scheme is powerful and robust 
enough to deal successfully with contaminants such as point 
sources, dust and the CMB signal. 
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5 CONCLUSIONS 



The filtering techniques described in this work have proven 
to be useful tools for the detection of discrete objects in 
multifrequency data, which are known to have a prescribed 
frequency dependence. The only other assumption about the 
objects is their spatial template; any other quantities of in- 
terest are directly calculated from the data. We have as- 
sumed in this work that the sources have circular symmetry, 
but the method can be extended to deal with asymmetric 
structures. Future work will take this point into account. 
We have applied these techniques to the problem posed by 
the presence of Sunyaev-Zel'dovich effect in multifrequency 
CMB maps. The best results are obtained with matched 
multifilters. The detection level of the method is compara- 
ble to other single-component separation methods, but the 
number of assumptions and previous data processing is re- 
duced to a minimum. Owing to the weakness of the signals 
that we aim to detect, the determination of the amplitude 
of the clusters is biased; this is the strongest limitation of 
the method. Therefore, the method should be complemented 
with a posteriori analysis (for example data fitting) in or- 
der to improve the amplitude estimation of weak clusters. 
Finally, the filters here presented can be used in other fields 
of Astronomy, such as X-ray observations, and large-scale 
structure detection. 
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